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We investigate the hydro-dynamic properties of a fluid simulated with a mesoscopic solvent model. 
Two distinct regimes are identified, the "particle regime" in which the dynamics is gas-like, and the 
"collective regime" where the dynamics is fluid-like. This behavior can be characterized by the 
Schmidt number, which measures the ratio between viscous and diffusive transport. Analytical 
expressions for the tracer diffusion coefficient, which have been derived on the basis of a molecular- 
chaos assumption, are found to describe the simulation data very well in the particle regime, but 
important deviations are found in the collective regime. These deviations are due to hydrodynamic 
correlations. The model is then extended in order to investigate self-diffusion in colloidal dispersions. 
We study first the transport properties of heavy point-like particles in the mesoscopic solvent, as 
a function of their mass and number density. Second, we introduce excluded-volume interactions 
among the colloidal particles and determine the dependence of the diffusion coefficient on the col- 
loidal volume fraction for different solvent mean-free paths. In the collective regime, the results are 
found to be in good agreement with previous theoretical predictions based on Stokes hydrodynamics 
and the Smoluchowski equation. 

PACS numbers: 47. ll.+j Computational methods in fluid dynamics 
82.20.Wt Computational modeling; simulation 
82.70.-y Disperse systems; complex fluids 



I. INTRODUCTION 

The dynamics of complex fluids such as colloidal sus- 
pensions, dilute or semi-dilute polymer solutions, biolog- 
ical macromolecules, membranes, and aqueous surfactant 
solutions, is often governed by the hydrodynamic behav- 
ior of the solvent. Due to a large separation of length 
and time scales between the atomic scale of the solvent 
molecules and the mesoscopic scale of the solute, direct 
simulation approaches with explicit atomistic solvent are 
prohibitively costly in computer time. Therefore, several 
mesoscale simulation techniques have been developed in 
recent years in order to bridge the length- and time-scale 
gap. In particular, lattice-gas automata (LGA) 0, 0, 
lattice Boltzmann methods (LBM) 0. [J ■ smoothed- 
particle hydrodynamics (SPH) 0, Q , dissipative particle 
dynamics mPD) [1H 

1 1 0| , direct simulation Monte Carlo 
(DSMC) [Hill, fluid particle dynamics [3, and oth- 
ers, have been investigated. The basic idea of all these 
approaches is very similar: To obtain hydrodynamic be- 
havior on length scales much larger than the atomic scale, 
the detailed interactions and dynamics of the solvent 
molecules are not important; instead mass and momen- 
tum conservation are the essential ingredients to obtain 
the correct hydrodynamic behavior. Therefore, the dy- 
namics on the microscopic scale can be strongly simpli- 
fied, as long as the conservation laws are strictly satisfied. 
The different methods listed above differ in the way the 
solvent dynamics is implemented. 

Two main classes of mesoscopic simulation techniques 
can be distinguished, which are lattice and off-lattice 
methods. Lattice gas and lattice Boltzmann methods fall 
into the first class, while direct simulation Monte Carlo, 
dissipative particle dynamics, and fluid particle dynamics 
fall into the second class. Off-lattice approaches have the 



advantage that Galilean invariance is typically satisfied. 
Moreover, the interaction of the off-lattice solvent with 
solutes such as colloids, polymers and membranes can be 
taken into account more naturally. 

The mesoscale simulation technique, which we are in- 
vestigating in this paper, was introduced by Malevanets 
and Kapral a few years ago. It is a variant of the 
DSMC method, in which binary collisions are replaced by 
multi-particle collisions in a prescribed collision volume. 
This method has been called multi-particle-collision dy- 
namics (MPCD) or stochastic rotation dynamics (SRD). 
It employs a discrete-time dynamics with continuous ve- 
locities and local multi-particle collisions. Mass and mo- 
mentum are conserved quantities and it has been demon- 
strated that the hydrodynamic equations are satisfied 
dill. 

Certain transport coefficients, in particular the viscos- 
ity, of this solvent model have been studied intensively. 
Analytical expressions have been derived from kinetic 
theory by generalizi ng p oint-like collisions to finite colli- 
sion volumes [iH ITtI Il8l Il9l | . The theoretical expressions 
describe numerical results very well. 

In this article, we study the transport coefficients as a 
function of the parameters of the MPCD fluid, in partic- 
ular the mean free path in units of the size of the collision 
volume. We find two distinct regimes, in which the dy- 
namics is either gas-like or fluid-like. This behavior can 
be characterized by the Schmidt number, which measures 
the ratio between viscous and diffusive transport. We 
find that MPCD allows us to tune the fluid behavior such 
that large Schmidt numbers are obtained and momentum 
transport dominates over mass transport. Analytical ex- 
pressions [171 Il8l Il9| for the tracer diffusion coefficient, 
which have been derived on the basis of a molecular- 
chaos assumption, are found to describe the simulation 
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data very well for large mean free paths, but fail in the 
fluid regime. The reason is a build-up of correlations 
among the fluid particles by hydrodynamic interactions, 
which leads to enhanced diffusion coefficients. We will 
show that the latter leads to non-exponentially decay- 
ing velocity-autocorrelation functions at small mean free 
paths. Independent of the mean free path, we find that 
the algorithm reproduces the algebraic long-time decay 
typical in fluids. 

In a further step, we investigate the diffusion of a heavy 
tracer particle in a MPCD solvent. It is very important to 
understand the contribution of the solvent dynamics on 
the solute diffusion. Two limiting situations are found: 
either Brownian or hydrodynamic behavior, depending 
on the collision time and the rotation angle. We explore 
the range of parameters where these different dynamical 
behaviors appear, and show how they emerge from the 
mesoscopic dynamics. 

Finally, we study self-diffusion in colloidal dispersions 
with excluded-volume interactions as a function of the 
volume fraction. To this end, the MPCD method is com- 
bined with molecular dynamic simulations. We find that 
such a hybrid model displays the proper dynamics for the 
same parameter regime where the hydrodynamic behav- 
ior is found for the fluid. Our results in the collective 
regime are in good agreement with previous theoreti- 
cal predictions based on Stokes hydrodynamics and the 
Smoluchowski equation 20]. 

II. THE MODEL 

The fluid is modeled by N point particles, which are 
determined by their positions rj and velocities Vj, with 
i = 1,...,N. Positions and velocities are continuous 
variables, which evolve in discrete increments of time. 
The mass m associated with the particles is taken to be 
the same, but more generally, different masses can be 
assigned. The algorithm consists of two steps, stream- 
ing and collision. In the streaming step the particles 
move ballistically according to their velocities during a 
time increment h, to which we will refer as collision time. 
Thereby, the evolution rule is 

n(t + h) = Ti(t) + hvi(t). (I) 

In the collision step, the particles are sorted into collision 
boxes, and interact with all other particles in the same 
collision box. The collision boxes are typically the unit 
cells of a d-dimensional cubic lattice with lattice constant 
a, although other geometries would be possible. The col- 
lision is then defined as a rotation of the velocities of all 
particles in a box in a co-moving frame with its center 
of mass. Thus, the velocity of the i-th particle after the 
collision is 

Vi(i + h) = v cm<i (t) + K(a) [ Vi (i) - v crrM (t)] , (2) 
where lZ(a) is a stochastic rotation matrix, and 
v C m,i(*) = J2'j' t \ mv j)/J2j m i s the velocity of the cen- 



ter of mass of all particles j, which are located in the 
collision box of particle i at time t. The conservation 
of local momentum and kinetic energy is guaranteed by 
construction. In two dimensions, the rotation of the rel- 
ative velocity is simply given by an angle ±a. Here a is 
a parameter of the model; the sign is chosen randomly 
for each cell. In three dimensions, various schemes for 
the random collisions are possible [lj, Il8l |2l| . The one 
employed in this paper consist in choosing a random di- 
rection in space for each box around which the relative 
velocities are rotated by an angle a. A detailed explana- 
tion of the implementation is given in Ref. |2l| . 

In order to ensure Galilean invariance for the full range 
of parameters, a random shift of the collision grid has 
to be performed in the execution of the collision step 
[lti l22j|. As a consequence of such a shift, the collision 
environment of each particle is independent of the aver- 
age local velocity, and no special reference frame exists. 
Random shifts also facilitate the transfer of momentum 
between neighboring particles. 

In the simulations, N particles are initially placed at 
random in a cubic system of linear extension L. The aver- 
age number of particles in a collision box is p = N(a/L) d , 
the scaled number density. Starting from an arbitrary 
distribution of velocities, only a few steps are required to 
reach the Maxwell Boltzmann velocity distribution. The 
equilibrium temperature T is then given by the average 
kinetic energy m (v?) = 3k B T, where ks is the Boltz- 
mann constant. In the simulations, we scale length and 
time according to x = x/a and t = t^/ksT /ma 2 , which 
corresponds to the choice m = 1, a = 1, and kgT = 1 
of reference units. The scaled mean free path is then 
given by A = h. Basic parameters and the definitions of 
dimcnsionless quantities are collected in Table [I] 



Parameters 
a Collision box size 
m Mass of the fluid particle 
T Temperature 
h Collision time 
a Rotation angle 
L Linear system size 
N Total number of particles 
g Mass density, g = Nm/L d 

A Mean free path, A = hsJksT /m 

DlMENSIONLESS QUANTITIES 

7 Decorrelation factor, 7 = (2/3)(l — cosa)(p — l)/p 
p Particles per cell, p = ga d /m = N(a/L) 

X Scaled mean free path, X = A/a 

Table I: Summary of relevant parameters for the simple fluid 
with the MPCD model. 



III. DYNAMICAL PROPERTIES 

The kinematic viscosity v = njp has been calculated 
theoretically O EE EE E3, EE Ea E3, H3| by means of 
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kinetic theory and its validity has been checked with sim- 
ulations. The total kinematic viscosity, v = u^ n + f C olh 
is the sum of two contributions, the kinetic viscosity Vkin 
and the collisional viscosity v co u , which have been calcu- 
lated in two and three dimensions. In three dimensions, 
the expressions 
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have been derived. 

The total kinematic viscosity has been determined nu- 
merically by the procedure explained in Ref. [24| • Briefly, 
a three-dimensional system is considered with periodic 
boundary conditions in two dimensions and planar walls 
in the third dimension. Stick boundary conditions at the 
walls are implemented by considering bounce-back colli- 
sions with the walls. A gravitational field is applied in one 
direction parallel to the walls. After a relaxation time, 
the system reaches a stationary state with a parabolic 
velocity profile between the walls and in the direction of 
the force. This is Poiseuille flow. It is known j2j| that the 
measured maximum velocity of the parabola is inversely 
proportional to the viscosity of the fluid. The viscosity 
data obtained in this way are presented in Fig.^together 
with the theoretical predictions of Eq. © . The obtained 
agreement is quite remarkable, in contrast to the case of 
other mesoscopic simulation techniques such as dissipa- 
tive particle dynamics 26 1. Density fluctuations can also 
be included in the theory [n| , which noticeably improves 
the agreement with the simulations results for small num- 
ber densities; for p = 5 and p — 10, these contributions 
are negligible. 

Alternative methods to determine the viscosity from 
simulations have been employed in Refs. 0] and [l7j . 
where a system under shear flow and vorticity correla- 
tions have been used, respectively. 

The ratio between the kinetic and the collisional con- 
tributions to the kinematic viscosity varies considerably 
with the model parameters, as can be seen easily from 
the theoretical expressions ®. In Fig. the total kine- 
matic viscosity and its two contributions are plotted as 
a function of the rotation angle and the collision time 
step. The collisional contribution is dominant for large 
collision angles and small collision times, while the ki- 
netic viscosity dominates in the opposite case of small 
collision angles and large collision times. 

Kinetic transport is due to the movement of the par- 
ticles themselves, i.e., when a particle moves it carries a 
certain amount of the relevant quantities as momentum 
and energy, while collisional transport is due to transfer 
of energy and momentum from one particle to another 
during collisions. In MPCD, kinetic transport is there- 
fore dominant when the mean free path is larger than the 
size of the collision box and for small values of the rota- 
tion angle. If the rotation angle is small, there is little ex- 




Figure 1: Dimensionless kinematic viscosity for the simple 
fluid in MPCD. Symbols are the simulation results, solid line 
is the total theoretical prediction, dotted line is the collisional 
contribution and dashed line the kinetic contribution. In both 
cases the system size is L/a = 20. In (a) a dependence is dis- 
played with A = 0.2 and p = 10. (b) shows the A dependence 
with a — 130 and p = 5. 



change of momentum between particles due to collisions. 
The situation where the kinetic transport dominates is 
characteristic for gases. In fluids the usual situation is 
the opposite, the transport of momentum is mainly due 
to collisions. 

A convenient measure of the importance of hydrody- 
namics is the Schmidt number Sc — v/D, where v is 
the kinematic viscosity and D the diffusion coefficient. 
Thus, Sc is the ratio between momentum transport and 
mass transport. It is known that this number for gases is 
smaller than but on the order of unity, while in fluids like 
water it is on the order of 10 2 to 10 3 . A prediction for 
the Schmidt number of a MPCD fluid can be obtained 
from the theoretical expressions for the kinematic vis- 
cosity, and the diffusion coefficient, see Eq. (|17(l below. 
In Fig. we plot the theoretical prediction for Sc as a 
function of the collision time for different values of the 
rotation angle. This shows that Sc becomes considerably 
larger than unity for the same range of parameters where 
the collisional viscosity is considerably larger than the ki- 
netic viscosity (Fig.^J. We will show that the dynamical 
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Figure 2: Theoretical Schmidt number versus collision time. 
The a and p parameters are specified in the plot. 



behavior in the two limits is fundamentally different. We 
will call the parameter region of large rotation angles 
and small collision times the collective regime, and the 
opposite region the particle regime. This classification 
has similar consequences as the one introduced in dissi- 
pative particle dynamics (DPD) |27j . although we do not 
investigate wave-length dependent properties here. 



IV. SIMPLE FLUID CORRELATIONS 

Correlations between particles are responsible for hy- 
drodynamic interactions. Therefore, we are interested in 
characterizing the velocity correlations in a MPCD fluid. 

A. Velocity Autocorrelation Functions 

An analytical expression for the velocity autocorrela- 
tion function (VACF) has been derived in Refs. [T^.[T^ |. 
The collision step in Eq. J2J can be rewritten as 

Vi(nh) =Vj((n - l)h)+ 

(K(a) - I) [vi((n - l)h) - v cra ,((n - l)h)] , 

(4) 

where I is the unit matrix and t = nh is the discretized 
time, with n the number of collision steps. By multi- 
plying this expression with the velocity at time zero and 
taking thermal averages, we obtain 



(Vi(nh)vi(0)) = (1 - 7tt ) (v t ((n - l)ft)Vi(O)) 

i 7a \ v cm, 

,((n- 1)^(0)}. 



(5) 



Here, the rotational average over an arbitrary vector A in 
three dimensions is obtained from geometrical arguments 
to be 



((K(a) - I) A) = --(1 - cosa) (A) = - 7q (A) . (6) 



This particular value of 7 Q arises from the implementa- 
tion of the rotation chosen in this paper. The remaining 
problem is to calculate the last term in Eq. J^J. First, we 
neglect density fluctuations in the average of the center of 

mass velocity, which yields (v cmj i(nh)) ~ (j^f^ v i^) I P- 
Furthermore, a molecular-chaos assumption implies that 



(vcm,i((n - l)/i)vi(0)) ~ - (v<((n - l)A)vi(O)) . (7) 

P 

This approximation means that of all the particles in the 
collision box of particle i after (n — 1) collisions, only 
particle i itself makes a non-zero contribution to the cor- 
relation function. This is the same as assuming that none 
of the other particles has any information about the state 
of particle i at any time. The correlation at a certain time 
step can then be expressed in terms of the previous time 
step as 

( Vi (n/i) Vi (0)) ~ (l - ^p7a) (vi((n- l)/i)vi(0)) ■ 

(8) 

This implies that in this approximation, the VACF shows 
an exponential decay, 

where the normalization factor follows from the equipar- 
tition theorem, (v?(0)) = 3k B T/m. The decorrelation 
factor 7 is defined as 
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^From Eq. @, a characteristic time tq = — /i/ln(l — 7) 
can be extracted. Up to this time, the VACF follows the 
exponential decay for every set of parameters. However, 
the collective phenomena responsible for the hydrody- 
namic behavior appear at much later times. 

In Fig. |3| simulation results of the VACF are presented 
for two different mean free paths A. The theoretical pre- 
diction 10 is also displayed for both values of A. For 
A = 1 the exponential decay is followed with very good 
accuracy until the crossover to the long-time tail behav- 
ior occurs. For A = 0.1 the purely exponential decay 
is followed only in the first collision; for long times, a 
long-time-tail behavior is observed similarly as for A = 1 . 
What is different in this case is that after the first colli- 
sion the system enters an intermediate regime where the 
VACF decay is significantly slower than the one described 
by the molecular-chaos approximation but is not yet the 
algebraic tail. Note that for the investigated rotation an- 
gle of a = 130, the mean free path A = 1 corresponds 
to the particle regime, while A = 0.1 corresponds to the 
collective regime. 

It is interesting to note that for short times, the VACF 
decays monotonically only in the case that the correla- 
tion parameter 7 is smaller than unity. If 7 > 1, Eq. @ 
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Figure 3: Normalized velocity autocorrelation function as a 
function of the dimensionless time for mean free paths A — 1 
and A = 0.1. Dashed lines correspond to the exponential 
decays in Eq. @. In both cases the number density is p — 5, 
the rotation angle a = 130, and the system size L/a = 20. 



predicts that the VACF exhibits damped oscillations. We 
have checked that this oscillatory behavior is indeed ob- 
served in the simulations. However, the viscosity curves 
show no particular features when this happens (compare 
Fig. ^ where the VACF for p = 10 becomes oscillatory 
for a > 132). 



B. Long-Time Tails 



It is well known 28, 29, 30] that the long-time behavior 
of the VACF in d-dimcnsional fluids in thermal equilib- 
rium shows a universal behavior. This corresponds to a 
power-law tail, for which the explicit form can be calcu- 
lated from a mode-coupling theory as |29|. 



C v (t) 



1 



d-1 



dp J [4tt(L> + v)t] d / 2 '' 



(11) 



where v and D are the transport coefficients of the fluid. 

The results obtained for the long-time behavior of the 
VACF are consistent with the general prediction for fluids 
in thermal equilibrium in Eq. (| 1 II) . The algebraic power 
£-3/2 j s c i ear i v reproduced in our simulations as can be 
seen in Fig. 0] The value of the amplitude in Eq. lllfl is 
related to the kinematic viscosity v and the diffusion co- 
efficient D. Since both values are known for the MPCD 
fluid and discussed in this paper, quantitative compari- 
son can also be performed. We find that the value for 
A = 1 is exactly reproduced by our simulations within 
the accuracy of the results, while the amplitude obtained 
for A = 0.1 is about 10% smaller than the theoretic pre- 
diction. Ihle and Kroll ^(| obtain good agreement in a 
two-dimensional MPCD fluid with the expected t -1 be- 
havior over a comparable time window. 

The effect of finite system size can be seen in Fig. 0] 
for times t > 20, where the VACF crosses over from the 
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Figure 4: Time dependence of the normalized velocity au- 
tocorrelation function. The parameters are the same as in 
Fig- El The data are compared with long-time tail prediction 
£-1/3 amplitude predicted in Eq. H1H is (within the 

statistical error) exact for A = 1 and about 10% larger for 
A = 0.1. 



algebraic to a faster, exponential decay. This effect is 
similar to that observed for the time dependence of the 
temperature autocorrelation function for a random-solid 
dissipative-particle-dynamics system 3 lj . There, it can 
be proved that the correlations decay faster after a time, 
where hydrodynamic modes become relevant which are 
truncated by the system size. 



C. Importance of Many-Body Correlations 



In the previous section, an exponential decay of the 
VACF has been theoretically predicted. This behavior 
is a consequence of the approximation in Eq. (JJJ which 
neglects any correlation among the particles in the same 
collision box at all times. In order to improve Eq. J7J, 
we have to go beyond the molecular-chaos approxima- 
tion. This is a formidable task. We start the procedure 
by calculating the center-of-mass correlation average for 
the first collision and, consecutively, the second and so 
on. For n = 1 the approximation in Eq. J7J is exact, 
(v c m,i(0)vj(0)) = (vf (0)) I p. This is the reason why for 
the first time step, C v (h) agrees perfectly in all simula- 
tions. For n = 2 it reads 
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where Ci denotes the number of particles that are neigh- 
bors of particle i at both times t = h and t = 0. We use 
the term "neighbors" for particles within the same colli- 
sion box. The approximation in Eq. (|7J) is recovered for 
(i = 1. This is the case when only the actual particle is 
considered to be in both collision boxes. As we have seen 
above, this is not a good approximation in the collective 
regime. 

We denote the average number of remaining neighbors 
that one particle is revisiting after n collisions as ( n . This 
number could in principle be calculated analytically by 
probabilistic arguments, but in order to get a flavor of the 
improvement that such numbers produce in the theory, 
we determine C, n numerically in our simulations. As ex- 
pected, these numbers strongly depend on the system pa- 
rameters. A detailed study has not been performed, but 
we have observed that the number of remaining neighbors 
seems to be a universal function of the root-mean-square 
displacement of the tagged particle. 
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Figure 5: Number of remaining neighbors in a collision box 
after n collisions as a function of the root-mean-square dis- 
placement. Solid symbols correspond to p — 5, a — 130 with 
mean free path specified in the legend. Open symbols cor- 
respond to p = 10, a = 110 with A = 0.6 (a) and A = 1 
(□). 



The measured numbers £„ are presented in Fig. [5J 
as a function of the root-mean-square displacement 
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(r(i) — r(0)) z ^ = \/6-Di, where D is the diffusion co- 
efficient and t = nh. The diffusion coefficient is the one 
obtained from the analytical expression which will be de- 
duced in the next section (see Eq. j!7))). The data for 
different mean free paths seem to fall onto a single mas- 
ter curve with reasonable accuracy. When the numerical 
values of D (discussed in the next section) are used in- 
stead of the theoretical result, the data collapse becomes 
even more accurate. For the large mean free path A = 1, 
the first collision takes place when \^6Dt/a ~ 2, which 
implies that d ~ 1 is a good approximation. Note that 
in the representation chosen in Fig. [SI £ n ~ 1 corresponds 
to the abscissa. The same displacement for a small mean 
free path A = 0.1 takes place when the particle has been 
involved in 80 collisions on average. The first collision 
for A = 0.1 takes place when the average displacement is 
much smaller and many of the particles are still in the 
same collision box, which makes £i ~ 1 a bad approxi- 
mation. Indeed, we can infer from Fig. [SJ that £i ~ 2.1 
for A = 0.1 and p = 5, and £i ~ 3.5 for A = 0.1 and 
p=10. 

Following the same procedure as employed in Eq. i|12[l , 
the velocity correlation function can be calculated for 
n = 3, 



M (2/Ov<(0)) = 



<^ 2 (0)> 
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where <5^ 2 is determined by 



C.2 + 5( 2 



k 2 (o)> 



(i,2) (1,0) 
3 I k 



(13) 



(14) 



This is the number of neighbors of particle i at the two 
times t = 2h and t = together with the neighbors of the 
neighbors, or the result of ring collisions. Let us consider 
two particles i and k, which are in the same collision box 
at t = 2h but not at t = 0. If one, k, has been neighbor 
of a third particle j at t = h and this j was neighbor of 
i at t = 0, then this combination also contributes to the 
correlation function. To obtain a reasonable prediction 
for this number is obviously not trivial. Furthermore, this 
relation will become more interconnected and difficult to 
predict for further time steps. It can be checked that with 
the approximations £„ = 1 and 5C, n — 0, Eq. i|13|) reduces 
to Eq. (JTJ, and consequently the exponential decay in 
Eq. © is recovered. 

Now we come back to the correlation average in Eq. (JSJ 
which can be expanded with the help of Eq. — with- 
out any approximation — 



(v i (n^)v i (0)) = <«?(0)>(l-7) n 
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The predictions for short times can be improved com- 
pared to Eq. 10 by employing the results of Eqs. 1(121) and 
(113(1 on the right-hand side of Eq. 1)15(1 , but setting £„ ~ 1 
and 5( n ~ for n > 3 as before. The result is shown in 
Fig. We observe that the prediction for the second 
collision C v (2h) now agrees perfectly with the simulation 
data, which confirms our arguments. Nevertheless, the 
prediction for further steps is still only a small improve- 
ment compared to the exponential decay in Eq. 
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Figure 6: Time dependence of the normalized velocity auto- 
correlation function. The dashed line is the exponential decay 
in Eq. crosses (x) are the simulation results and pluses 
(+) are the predicted values obtained by employing the £ n 
numbers, as indicated in Eq. (1151 . The parameters of the 
simulation are p = 10, a = 110, A = 0.1 and L/a = 20. The 
inset is a zoom into the regime of the first few collisions. 

The most relevant conclusion at this point is that in 
the collective regime the MPCD algorithm accounts for 
many-body collisions which are crucial for the build-up 
of correlations. This is known to be the origin of the 
hydrodynamic behavior in fluids. 



V. SELF-DIFFUSION 

We study now the consequences of the different behav- 
ior in the two hydrodynamic regimes, which have been 
introduced in Sec. EI on the self-diffusion coefficient. 



A. Diffusion Coefficient 

In the Green-Kubo formalism, the self-diffusion coeffi- 
cient is given by D = | J Q dt (v(t)v(O)) . In the case that 
the time is discretized the integral has to be replaced by 

mn 



^=3 



^ 2 (o))/>- 



£ (y(#(0)) 



(16) 



In order to obtain an analytical prediction for the diffu- 
sion coefficient, an expression for (v(n/i)v(0)) is required. 



The Brownian approximation for the VACF given by 
Eq. © yields 



m \7 2 



(17) 



with 7 defined in Eq. {TDJ. This expression coincides with 
that of Ref. with a different notation. 
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Figure 7: Relative deviation AD = (D S i m — Do) /Do of the 
simulated diffusion coefficient from the Brownian approxima- 
tion, as a function of the scaled mean free path A. Full circles 
are simulation results, the solid line represents the analyti- 
cal expression of Do in Eq. 1)170 . Simulation parameters are 
a = 130, p = 5 and L/a = 20. 

In the simulations, the diffusion coefficient is deter- 
mined by a linear fit of the mean-square displacement 
for long times. We have checked that equivalent results 
for D are also obtained directly from the VACF by em- 
ploying Eq. 1 )16( 1 . 

Fig. [7] shows the relative deviation AD — (D S i m — 
Dq)/Dq of the diffusion coefficient from the expression 
(117)1 . This expression should be a good approximation 
as long as the exponential decay 10 of the VACF ap- 
plies. This is indeed the case for A > 0.6, what means 
that the long-time tail for these values has a negligible 
contribution for the diffusion coefficient. This is reason- 
able since the deviation from the exponential behavior 
appears when the VACF has decayed typically by three 
orders of magnitude (see A = 1.0 in Fig- - I n contrast, 
Fig. [7] shows that the deviation from the Brownian be- 
havior ()17(l increases with decreasing A for A < 0.5. This 
can be understood from the VACF since for small A the 
deviation from the exponential decay appears much ear- 
lier. Fig.|21shows that for A = 0.1 the VACF has decayed 
only by about one order of magnitude when the devia- 
tion starts. This translates into a noticeable increment 
of the diffusion coefficient. This difference can be under- 
stood as a hydrodynamic enhancement of the diffusion 
coefficient for large values of the Schmidt number. 

The diffusion coefficient for a simple MPCD fluid in 
two dimensions has been determined by Ihle and Kroll 
In their Fig. 15, results for A = 0.113 are pre- 
sented as a function of the rotation angle; deviations from 
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the theoretical prediction are found for large values of a, 
which is in the range of parameters which we identify as 
the collective regime. They arrive at a similar conclusion 
that this is due to multiple encounters among particles. 
In three dimensions, some numerical results of the dif- 
fusion coefficient have been presented in Ref. and 
good agreement with the molecular-chaos approximation 
has been found for a large range of number densities. 
However, the employed parameters (which correspond to 
A > 0.5) all belong to the particle regime, where we ar- 
gue that a good agreement with the theory should be 
expected. 

At this stage we come back to the discussion in Scc. lIIII 
about the Schmidt number Sc = v j D. The analytic ex- 
pression can be calculated from the viscosity v in Eq. @ 
and the diffusion coefficient in Eq. (|17|l . as was already 
pointed out in Refs. ^7ll3llp|. Note that Sc increases 
rapidly for small values A < 1 of the mean free path, 
where Sc ~ tiT 2 . This allows arbitrary large values of the 
Schmidt number. Although very small values of the colli- 
sion time significantly reduce the efficiency of the simula- 
tions, there is a range of A- values which are not too small 
but still display fluid behavior corresponding to high Sc. 
On the other hand, the hydrodynamic enhancement of 
the diffusion coefficient in the collective regime leads to 
values of Sc which are smaller than predicted by the an- 
alytical approximation. By substituting the numerically 
determined diffusion coefficient, it can be checked that 
Sc is indeed smaller, but still large enough to display a 
fluid-like behavior. 



B. Continuum Time Limit 

It is interesting to discuss the limit of small collision 
times h — > 0, and small rotation angles a — ► 0. The 
leading contributions in the theoretical expressions © 
of the kinetic and collisional viscosity read in this limit 



Vcoll 
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36a 
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« 3 7p 



with 7 P defined in Eq. (|1U|) . This result shows that a 
finite viscosity is obtained in the continuum limit only if 
the ratio a 2 /h is kept constant. The additive term due 
to discrete times in Eq. I|17|) naturally vanishes in the 
continuum limit, because 7 ~ a 2 . 

The expressions (|18|) for the kinetic and collisional con- 
tributions to the viscosity show that the collective regime, 
where v co u 3> Vkim corresponds to a 2 /h 3> 1 in the con- 
tinuum limit. In this regime, the leading contribution to 
the diffusion coefficient <|17ll is found to be 



D 



3k B T ( h 

7 P 



The related Schmidt number 

1 



2 2 
171 % 



108 ak B T 



(19) 



(20) 



can be very large since a 2 /h ^> 1. This shows that the 
model has a proper continuum limit. However, due to 
the requirement of very small collision times, this limit is 
not very convenient from a computational point of view. 

It is very satisfactory to see that the Stokes-Einstein 
relation is satisfied in this case, since the diffusion coeffi- 
cient is inversely proportional to the viscosity, 



D 



6npv cM R 



with 



(21) 



and defines an effective particle radius inversely propor- 
tional to the number density. We want to emphasize, 
however, that the Stokes-Einstein relation is not only 
satisfied in the continuum limit, but always when the 
additive term 1/2 in Eq. 117(1 can be neglected and the 
collisional dominates the kinetic viscosity. In this case, 
Eq. l|?T|) is also valid. 



VI. DYNAMICS OF EMBEDDED PARTICLES 

After the behavior of a simple MPCD fluid has been 
characterized, the next important question is how com- 
plex fluids can be modeled. As first step, we investi- 
gate the behavior of a single heavy point-like particle, 
which could represent a solute particle or a colloidal 
sphere embedded in a simple fluid. Also, the monomers 
in a polymer chain can be represented as point parti- 
cles |3J, |35j, |36|, |3jj • This is a quite convenient strategy, 
since the solute-solvent interactions are modeled by just 
including the point-like solute particles in the collision 
step. Then we study different concentrations of these 
heavy particles. 



A. Single Heavy Tracer Particle 

For the simulation of heavy point-like particles em- 
bedded in a solvent, the algorithm is the same as de- 
scribed for the simple fluid in Sec. [H] The only point 
where the higher mass plays a role is in the calcula- 
tion of the velocity of the center of mass, where the 
different particle masses have to be taken into account 



In thermal equi- 



via v cm) j(t) = J2j( m j v j)/J2 
librium, the average kinetic energy of light and heavy 
particles is the same. Therefore, the average momentum 
of the heavy particle of mass M is a factor (M/m) 1 / 2 
larger than the average momentum of a light particle. 
This implies that a heavy particle has a larger contribu- 
tion in the center-of-mass velocity than a light particle. 
Since the center-of-mass velocity and therefore also the 
velocities of all particles after the collision step depends 
on M and the mass mp of the solvent particles in a col- 
lision cell, the effective coupling between the solvent and 
the solute must depend in general on the ratio M/(mp). 

We denote the heavy particle position and velocity 
with capital letters R and V. Of course, all types of 
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particles are involved in the center-of-mass calculation 
or other sums over particles. The VACF can be calcu- 
lated in the molecular-chaos approximation as explained 
in Sec. IIV Al except for the center-of-mass correlation in 
Eq. J7J), which for the heavy particle yields 



(V cm ((n - l)h)V(0)) 



M 



mp 



M 



(V((n - l)fr)V(0)> 



(22) 

because in the collision box of the heavy particle the total 
mass is (M + mp). The correlation at time zero depends 
now on the heavy particle mass, 



<V») 



' m ' 



(23) 



By inserting these results in the expression equivalent to 
Eq. (jHJ , we obtain the molecular-chaos approximation for 
the normalized VACF of the heavy particle, 



_ (V(nfe)V(O)) ^ 

Cv(t) = (vm) ( 7) ' 

where the decorrelation factor 7 is now given by 

mp 



1 = lo 



mp + M 



la 71, 



(24) 



(25) 



and 71 is defined for one heavy particle in the presence 
of p fluid particles, in contrast to 7 P in Eq. (frU|) , where a 
fluid particle is surrounded by (p— 1) other fluid particles. 
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Figure 8: Time dependence of the normalized velocity auto- 
correlation function for different heavy particles. Simulation 
parameters are A = 0.1, a = 130, L/a = 20 and p = M/m. 
Dashed lines are simulation results and the solid line is the 
molecular-chaos approximation 1241 1. 

In Fig. [S] results for the normalized VACF of one heavy 
particle in the collective regime are presented for differ- 
ent values of its mass. The solvent mass density has been 
chosen to be equal to the solute mass, i.e., p — M/m. In 
this way, 7 = 7 Q /2 and the analytical expression l|24|l is 
independent of the heavy particle mass. Fig.[S]shows that 
after the second collision all the simulation data exhibit a 
non-exponential decay. This is not very surprising, since 



a similar behavior was observed for the simple fluid in 
Fig. |3 for parameter values within the collective regime. 
A slightly slower decay is displayed at lower number den- 
sity p, but an asymptotic curve is clearly approached for 
large values of p. The deviations for small p are due to 
the presence of density fluctuations. 
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Figure 9: Time dependence of the normalized velocity auto- 
correlation function of a heavy particle of mass M — 5m for 
mean free path A = 1 and A = 0.1. Dashed lines correspond 
to the exponential decays in Eq. (1241 . In both cases the num- 
ber density is p — 5 and the rotation angle a = 130. Compare 
with Fig. ||l 

The dependence of the VACF of a single heavy tracer 
particle of mass M = 5m on the mean free path A of 
the solvent is shown in Fig. |5J Corresponding results 
for the simple fluid are shown in Fig. [3] For A = 0.1, 
the qualitative behavior of tracer particles with M = m 
and M = 5m is very similar. The first collision perfectly 
follows the molecular-chaos approximation, followed by a 
slower-than-exponential decay for intermediate times and 
a crossover to a power-law decay for long times. However, 
note that since the exponential decay is slower for the 
heavy particle, the deviations from Brownian behavior 
appear when the VACF has decayed to approximately 
one third of its original value for the employed values of 
p and a, while for the simple fluid case the VACF has 
decayed to 6% of its original value. This implies that 
the hydrodynamic enhancement is more pronounced for 
particles of larger mass. For A = 1 , small deviations from 
the exponential decay are visible for short times; for long 
times, the crossover to the power-law behavior can be 
seen. 

Analytical approximation for the diffusion coefficient 
can be calculated similar to Sec. IV Al It reads, 



Dn 



k B T u fl 1 



(26) 



where the decorrelation factor 7 is now given by Eq. 12511 . 

Simulation results for the diffusion coefficient Dm of a 
heavy tracer particle are plotted in Fig. 1101 as a function 
of the mass M/m, for fixed solvent density p = 5 and 
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Figure 10: Relative deviation of the simulated diffusion coeffi- 
cient Dm from the Brownian approximation Do, in Eq. 1261 . 
as a function of the heavy particle mass for p — 5. This 
deviations represent the hydrodynamic contribution to the 
diffusion coefficient Dh = Dm — Do in units of the Brownian 
contribution. Symbols are simulation results and the dashed 
fine is a guide to the eye which represents a 75% enhancement 
of the hydrodynamic term over the Brownian one. 
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Figure 11: Hydrodynamic contribution to the diffusion co- 
efficient in units of the Brownian contribution as a function 
of the scaled mean free path A. Symbols are simulation mea- 
surements, and the ordinate zero axis represents perfect agree- 
ment with the analytical expression Do. Simulation parame- 
ters are a = 130, p = 5 and L/a = 20. Compare with Fig. |7| 



two different sets of parameters. The agreement of the 
simulations with the approximation i|26|) is again very 
good for parameter values within the particle regime, A = 
1 and a — 45, but not within the collective regime, A = 
0.1 and a = 130. This is the same behavior as observed 
in the simple fluid (see Fig. [7J| and indicates again the 
presence of a hydrodynamic contribution to the diffusion 
coefficient in the collective regime. 

In Fig. ^2 the hydrodynamic contribution to the dif- 
fusion coefficient (in units of the Brownian contribution) 
is plotted as a function of the scaled mean free path A for 
a heavy tracer particle of mass M = 5m and for a sim- 
ple fluid tracer particle (compare Fig. El ■ It can be seen 



that Dh increases considerably for small A in both cases. 
This increment is significantly more pronounced for the 
heavy particle, which corresponds to the slower decay of 
the VACF in Fig-Ofor the larger mass. A small deviation 
of the VACF from the exponential decay was observed in 
Fig. at short times for A = 1. This deviation translates 
into the small hydrodynamic enhancement of the diffu- 
sion coefficient of the heavy particle that can be seen in 
Fig. 1111 even at "large" mean free paths A ~ 1. 

Fig.llOlshows that for a fixed density p in the collective 
regime, the hydrodynamic enhancement increases with 
increasing mass of the solute particle until M/m ~ 2p, 
and then levels off and becomes independent of the solute 
mass for M/m ^> p. This is consistent with the diffusion 
behavior of colloidal spheres, where the diffusion coeffi- 
cient is independent of the mass of the colloidal particles. 

Kikuchi et al. [l9| determine numerically the friction 
coefficient acting on a particle of mass M and velocity v 
in a MPCD solvent. Their simulation results, for a fluid 
of A ~ 0.9, compare nicely with the analytical prediction, 
independently on the mass of the particle. However, we 
want to point out that this agreement is not very surpris- 
ing, since their result is obtained from the velocity au- 
tocorrelation function after the first collision step, where 
the molecular-chaos approximation is always exact (see 
Sec. ESt - 

The increase of the hydrodynamic coupling of solute 
and solvent with increasing solute mass can be under- 
stood as follows. The relative mass of the solute and 
solvent particles appears in the collision step via the cal- 
culation of the center-of-mass velocity. If solute particles 
have the same mass as solvent particles and there is a 
large number of solvent particles per cell, the solvent par- 
ticles transfer a large random momentum to the solute 
particle. Simultaneously, the effect of the solute particle 
momentum on the solvent is small. For this reason, the 
hydrodynamic contribution to the diffusion constant of a 
particle of equal mass, shown in Fig. is only of the or- 
der of 30% for the largest Schmidt number considered. In 
contrast, this hydrodynamic enhancement is 65% when 
M/m ~ p and 75% when M/m > 2p, as can be seen in 
Fig. A very large mass of the solute particle is not 
very convenient either, because it implies a large ballistic 
regime and a long diffusion time. Therefore, we conclude 
that a mass M/m ~ p for the solute particle is a optimal 
choice to enhance the hydrodynamic coupling between 
solute and fluid particles. 



Finite Concentration of Heavy Point-Like 
Particles 



At a finite concentration of solute particles, an im- 
portant question is to which extent solute particles build 
up hydrodynamic interactions among themselves through 
the fluid particles when simulated with MPCD. We study 
therefore systems with different concentrations of heavy 
particles for sets of parameters within the particle and the 
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collective regimes, respectively. We address this question 
by investigating the tracer-diffusion coefficient. 

Simulations with different heavy particle concentra- 
tions are performed by changing the total number Nm 
of heavy particles but keeping fixed the volume V = L 3 
and the number N of solvent particles. The corre- 
sponding number density of heavy particles is defined 
as <f> = Nm{o,/L) 3 . In Fig. the diffusion coefficients 
for three different values of the mean free path are dis- 
played. Very surprisingly, when the data are normalized 
by the corresponding diffusion coefficients in the limit of 
vanishing density (j), all three data sets, which are both in 
the particle and the collective regime, collapse onto a sin- 
gle curve. We recall that the hydrodynamic enhancement 
for the diffusion coefficient of a single heavy particle, here 
denoted as -Dm(O), is quite different among these three 
values of A (see Fig. Illf) . It can be inferred from the 
data collapse in Fig. ^] that there is no extra hydrody- 
namic contribution among these heavy particles, which 
is consistent with the idea that there is no hydrodynamic 
screening for point particles [H, |3|j . 
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Figure 12: Diffusion coefficients for a heavy particle as a func- 
tion of the concentration <f> = Nm((i/L) 3 , normalized with the 
diffusion coefficient Dm (0) of heavy particles at zero num- 
ber density. The dashed line is the analytical approximation 
from Eqs. 12611 and 1291 . symbols correspond to the simulation 
data with A specified in the legend. The other parameters are 
M/m = p = 5, a = 130 and L/a = 20. 

The dependence of the diffusion coefficient on the 
heavy particle number density can be understood along 
the same lines as for the simple fluid or the single heavy 
particle. We assume that in each collision box there is 
a fixed number of fluid particles p, but that the number 
of heavy particles n, fluctuates from one collision box to 
another. The probability P(n) of a given heavy particle 
to be found in a cell with a total of n — 1 other heavy 
particles is given by the Poisson distribution function, 
P(n) = e~"^ n_1 /(n — 1)!. The corresponding decorre- 
lation factor for a heavy particle in a collision box with 
(n — 1) other heavy particles and p fluid ones is 



compare the definition of 71 in Eq. I|25[l for a single heavy 
particle in a collision box. The diffusion coefficient is 
then given by Eq. (|26|) . where the decorrelation factor is 
now 7 = 7 Q Y] P(n)"f n . In the regime of low number 
density, <f> <C 1, this implies 



7 = 7« [(1 - 4>)ll + 072 + O (4> 2 )} . 



(28) 



In the special case of p — M / m, the sum can be evaluated 
analytically and yields 



7 = 7c (l-(e-* + ^-l) /0 2 ). 



(29) 



In Fig.^]the simulation data for the normalized diffu- 
sion coefficient at different volume fractions are compared 
with the theoretical prediction obtained from Eq. 12(j|) 
with the decorrelation function in Eq. (|29|1 . It can be 
seen that this prediction overestimates the values for the 
diffusion coefficients. Further studies are required to un- 
derstand the origin of this deviation. 



VII. HYBRID DYNAMICS 

In order to go one step further in the development of 
an efficient simulation technique for suspensions of col- 
loidal particles with MPCD, we next investigate the ef- 
fect of excluded-volume interactions between the heavy 
particles. To this end, the MPCD algorithm has to be 
combined with standard molecular dynamics (MD) for 
the solute particles. 



A. The Model 

We consider a dispersion of spherical colloidal parti- 
cles in three dimensions. The interactions of solvent par- 
ticles among themselves and with colloids take place in 
the MPCD collisional step, exactly in the same way as 
described for the heavy point-like particles in Sec. IV! Al 
However, the streaming step Q is used only for the sol- 
vent particles. The position update of the colloidal par- 
ticles is performed in several MD steps between MPCD 
collisions. In these MD steps, colloids interact via an 
excluded-volume potential. We use the truncated repul- 
sive Lennard Jones potential |4Cj 



V 



RLJ 



(r) 



4e 



■e, r < r n 
0, r > r, n 



(30) 



7n = 1 - M/(pm + nM), 



(27) 



where r is the distance between the centers of the col- 
loidal particles. The parameter a is related to the parti- 
cle diameter; it is chosen to equal the collision box length, 
a = a, so that there is typically no more than one colloid 
particle in each collision box. The potential strength is 
taken to be equal to the thermal energy, e = kgT, the 
cut-off radius is r m i n = 2 1 / 6 <r, and the mass of the par- 
ticles is taken to be M = 5m. The MD time steps are 
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integrated with the velocity- Verlet algorithm 4d| with a 
time step At = 0.002,/e/ma 2 . 

In other words, we consider a system of colloidal parti- 
cles interacting through repulsive Lennard Jones poten- 
tials whose positions and velocities evolve in discrete time 
intervals At. This procedure is interrupted every h/ At 
steps for the interaction with the fluid particles. This 
interaction is a MPCD event where solvent and solute 
particles interchange momentum. This implies that the 
solvent particles can enter the cores of the colloidal par- 
ticles, but the colloids cannot interpenetrate each other. 

The hybrid model described here is a variant of 
the model introduced previously by Malevanets and 
Kapral [TEl 142^ . In their model, both the solute-solute 
and solute-solvent interactions were taken into account 
through excluded- volume potentials with MD, and only 
the solvent-solvent interactions were mesoscopically de- 
scribed through MPCD. The advantage of the model de- 
scribed here comes from the fact that in the MD steps 
just the solute particles are considered. This leads to a 
considerable speed up of the simulations. 
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Figure 13: Diffusion coefficient for colloidal particles without 
solvent as a function of the volume fraction (p. Symbols are 
simulation results, the dashed line corresponds to the ana- 
lytical prediction 132H . The inset is a zoom over the small 
values of the diffusion coefficient, and the solid line is a linear 
extrapolation for large values of ip. 



B. Diffusion in Colloidal Dispersions 

We measure the diffusion coefficient of the dispersion 
through the mean-square displacement of a tracer par- 
ticle, as before. Simulations are performed for different 
colloidal concentrations. The volume fraction of colloidal 
particles (p is the fraction of the total volume V occu- 
pied by the colloidal particles, ip = (tt j%)a\^pMi where 
the effective diameter <J e j f is determined by the Barker- 
Henderson expression |43( 

<?eff = f dr [1 - exp(~V RLJ (r)/k B T)] . (31) 
Jo 

For our choice of Lennard- Jones parameters, this gives 
(Jeff — l.Olcr. The number density of colloidal particles 
is pm = (Nm ~ — Nm/V, where Nm^s the number 
of heavy particles with excluded- volume interactions. 

For later comparison and better understanding of our 
hybrid model results, we recall first the basic behavior 
of a system with excluded-volume interactions only. In 
Fig-EH we show the results for the diffusion coefficient of a 
MD simulation of repulsive Lennard- Jones particles. Ki- 
netic theory for hard spheres predicts in the low-density 
limit H3 

Dmd(^) = A-J^—, (32) 
8cr^ // V ttM p M 

This analytical prediction is depicted in Fig. 1131 together 
with the simulation results. It can be seen that for 
small volume fraction, the ip^ 1 behavior is properly re- 
produced, while for large volume fractions a linear be- 
havior can be inferred. 

The density dependence of the self-diffusion coefficient 
of colloidal hard spheres in a hydrodynamic bath has 



been calculated in Ref. pfj|. 

Ds(p) = D s (0) [l - 2.1ip + O (ip 2 )] . (33) 

The diffusion coefficient now decreases linearly with the 
volume fraction, in contrast with the kinetic theory re- 
sult for a gas of hard spheres. In the calculation of 
Eq. 1|33|) . Brownian and hydrodynamic terms have to be 
considered, and it has been found that the hydrodynamic 
terms almost cancel. For a colloidal dispersion in a Brow- 
nian bath |2fJ the first-order correction in Eq. I|33|l equals 
— 2.0ip. Thus, no significant differences are expected be- 
tween Brownian and hydrodynamic measurements of the 
diffusion coefficient. 

Simulation results with the hybrid method are shown 
in Fig. 1141 The simulations presented here are performed 
with rotation angle a — 130, fluid number density p = 5, 
and mass M — 5m of the colloidal particle. We vary the 
mean free path between A = 0.02 and A = 2.0. 

In the limit of very small volume fractions, the repul- 
sive interactions between colloids are negligible, and the 
colloidal dispersion will behave as the dispersion of heavy 
point-like particles presented in Sees. IVlXl and IVIBI In 
this limit, we know from Eq. (|26|l that the diffusion co- 
efficient -D(O) increases with the mean free path A (with 
-D(O) ~ A in the molecular-chaos approximation). The 
decrease of D(ip) with decreasing A displayed in Fig. 1141 
arises then as a natural consequence. Furthermore, the 
MPCD interactions of the colloid particles with the fluid 
imply that the self-diffusion coefficient at small densi- 
ties does not diverge but goes to finite value dictated by 
Eq. d§. 

For small but finite volume fractions, in the case of 
large values of A, we observe a behavior reminiscent of 
the ip -1 decay of hard-sphere gases, instead of the linear 
decrease expected from Eq. (|33f) . This can be understood 
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Figure 14: Diffusion coefficient as a function of the volume 
fraction of colloidal particles dispersions interacting with a 
solvent represented with MPCD at different collision times. 
For comparison also plotted the MD results of Fig. [T31 



since, in the limit of very large mean free paths, the col- 
loids will essentially interact with each other rather than 
with the solvent. This behavior is not seen in experiments 
of colloidal dispersions, because the diffusive length scale 
is typically much smaller than the diameter of the parti- 
cles. 
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Figure 15: Dependence of the normalized diffusion coefficient 
on the volume fraction tp of colloidal particles. The same 
data are shown as in Fig. 1141 The normalization factor D(0) is 
obtained by extrapolation of the data to zero volume fraction. 
The solid line corresponds to the hydrodynamic prediction in 
Eq. 

Therefore, the appropriate parameters for the model- 
ing of colloidal dispersions have again to be chosen in the 
collective regime. In Fig. 1151 the normalized diffusion co- 
efficient is shown, where D(0) is extrapolated from the 
simulated data. The linear behavior in Eq. (|33l) is in- 
deed observed for the smallest values of the mean free 
path, A = 0.02 and A = 0.1, within the accuracy of the 
simulations. Thus, we find that in order to obtain the 
theoretically predicted behavior l|33|) from simulation of 
the MD - MPCD hybrid model, small values of the mean 



free path and large values of the rotation angle a are 
required, i.e. parameters in the collective regime. 

However, an almost identical dependence of the diffu- 
sion coefficient on the volume fraction is predicted theo- 
retically in the absence of hydrodynamics interactions. In 
order to investigate this point in more detail, we have per- 
formed simulations of a hybrid model similar to the one 
presented here, but with a completely Brownian solvent. 
One way of transforming a MPCD fluid in a Brownian 
solvent has been introduced by Kikuchi et al. , where 
the velocities among all the fluid particles are randomly 
interchanged after each MPCD collision step. We pro- 
pose an alternative method which does not consider any 
solvent particles. Instead, at every h/At steps the MD 
dynamics is interrupted for a rotation of the (full) veloc- 
ity of each colloid around a random axis by and angle a. 
In this case, the diffusion coefficient at zero volume frac- 
tion -D(O) is given by Eq. I|26|) but the decorrelation factor 
is 7 = 7 Q with 7„ of Eq. ©. The simulation results of 
the dependence of the diffusion coefficient on the volume 
fraction are quite similar to those displayed in Fig. ED 
The data for D((p) follow a linear decay only for very 
small values of A, where the friction is large and D(0) is 
small enough to represent a fluid. For large values of A, 
D(ip) has a concave shape, reminiscent of the </? _1 behav- 
ior of gases, similarly as observed for the hydrodynamic 
simulations. 

Simulations with a similar hybrid method of a two- 
dimensional colloidal suspension have been reported by 
Falck et al. [3j| . In the majority of the presented re- 
sults, they consider excluded-volume interaction among 
colloids but not between colloids and solvent particles. 
They measure an apparent tracer diffusion coefficient 
(since the diffusion coefficient in two dimensions diverges 
with increasing system size) of the colloids for differ- 
ent concentrations. Three different Schmidt numbers are 
studied. A similar trend in the data is observed as in our 
simulations: the normalized diffusion coefficient increases 
with increasing Schmidt number, in particular for inter- 
mediate values of the volume fraction (p. However, our 
interpretation is different. While Falck et al. attribute 
this effect to hydrodynamics, we believe that it is due 
to the crossover from gas-like to diffusive behavior of the 
colloidal dynamics. 

In summary, our hybrid model describes the dynam- 
ics of a dispersion of hard-sphere colloids very well in 
the collective regime of the solvent. In the hydrodynamic 
interaction, only the leading contribution for large dis- 
tances is included in our model. This implies that lu- 
brication forces between neighboring particles at short 
distances, as well as the coupling between rotational de- 
grees of freedom, are neglected. We conclude from the 
very weak dependence of our results for the normalized 
diffusion coefficients on the mean free path, which con- 
trols the strength of the hydrodynamic interaction, that 
our model works very well for not too concentrated col- 
loidal dispersions. 
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VIII. SUMMARY AND CONCLUSIONS 

In this paper we have performed a detailed analysis 
of the hydrodynamic properties of a fluid simulated with 
MPCD. We identify two hydrodynamic regimes in terms 
of the parameters of the MPCD algorithm. The parti- 
cle regime is characterized by dynamical properties be- 
ing closer to those of a gas than to those of a liquid. 
The Schmidt number is small and the dominant transport 
mechanism is kinetic transport. This is the regime ob- 
tained for large values of the collision time and/or small 
values of the rotation angle. The second and more rele- 
vant regime for fluid simulations is the collective regime. 
In this regime the Schmidt number is large and colli- 
sional transport dominates over kinetic transport — this 
characterizes liquid-like behavior. These properties are 
obtained for large values of the rotation angle and small 
values of the collision time. 

Different quantities have been measured in both 
regimes. The main conclusion is that the diffusion co- 
efficient shows a hydrodynamic enhancement in the col- 
lective regime. In the study of the VACF we observe that 
the behavior can be understood in both regimes as an ex- 
ponentially decay for short times and algebraic decay for 
long times. In the particle regime, a simple crossover be- 
tween both behaviors is observed while an extra interme- 
diate behavior is displayed in the collective regime. This 
intermediate behavior of the VACF is typically a slower 
than the initial exponential decay. We have shown that 
the origin of this intermediate decay region is due to the 
build-up of correlations by many-body collisions, which 
is in conceptual agreement with the hydrodynamic be- 
havior. The theoretical predictions for the diffusion coef- 
ficient are based on a molecular-chaos assumption, which 
gives an exponential decay of the VACF. Consequently, a 
deviation from the theoretical prediction is found in the 
collective regime. This deviation can be understood as a 
hydrodynamic contribution to the Brownian value. 

In a further step, we have investigated the differences 
between the particle and the collective regime for com- 



plex fluids. We have studied the behavior of heavy parti- 
cles embedded in the MPCD fluid which can represent so- 
lute or colloidal particles dissolved in a simple fluid. This 
study demonstrates that optimal hydrodynamic coupling 
occurs when the mass of the tagged particle is on the or- 
der of the solvent mass in a collision cell. 

In order to describe colloidal dispersions at finite vol- 
ume fractions, it is necessary to account for excluded vol- 
ume interactions among colloidal particles. To this end, 
a hybrid model was studied, which combines MPCD for 
the solvent with MD simulations for the colloidal parti- 
cles. We show that only for parameters within the collec- 
tive regime does the hybrid model reproduce the proper 
hydrodynamic behavior. In this case, the results agree 
well with the theoretical calculations with and without 
hydrodynamic interactions, as well as with experimental 
results. 

A more precise modeling of colloidal particles would 
require new interactions among fluid and colloids such 
that fluid particles would not freely travel through col- 
loidal particles and eventually angular momentum could 
be interchanged among them. 

In the future, it will be interesting to explore in which 
applications a more detailed description of colloidal inter- 
actions is necessary, as compared to our simplified model 
which allows more particles and larger system sizes, and 
is therefore well suited to study cooperative phenomena. 
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